#include "include/MultigridSolver.h"
#include "include/EquationSolver.h"

using namespace std;

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-31 12:37
 *
 * @description: 测试泊松方程
 */

#define Dim 2
#define PI 3.14159265358979

Vector u(Vector x)
{
    return sin(2 * PI * x[0]) * sin(2 * PI * x[1]);
}

Vector du1(Vector x)
{
    return 2 * PI * cos(2 * PI * x[0]) * sin(2 * PI * x[1]);
}
Vector du2(Vector x)
{
    return 2 * PI * sin(2 * PI * x[0]) * cos(2 * PI * x[1]);
}

Vector f(Vector x)
{
    return -4 * PI * PI * (sin(2 * PI * x[0]) * sin(2 * PI * x[1]) + sin(2 * PI * x[0]) * sin(2 * PI * x[1]));
}

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-31 13:54
 *
 * @description: 狄利克雷边值/周期边界
 */

Real testPoisson(int n)
{
    Domain<Dim> domain({0, 0}, 1);

    // 边界条件
    BoundaryCondition<Dim> bc(f, domain);
    bc.alpha_[0] = 1;
    bc.alpha_[1] = 0;
    bc.alpha_[2] = 1;
    bc.alpha_[3] = 1;
    bc.beta_[0] = 0;
    bc.beta_[1] = 1;
    bc.beta_[2] = 1;
    bc.beta_[3] = 0;

    // 网格
    Grid<Dim> g(n + 1, 0);

    // 导函数
    FuncX Du[Dim];
    Du[0] = du1;
    Du[1] = du2;

    // 填充边值函数
    bc.fillPoisson(u, Du);

    // 多重网格求解器
    MultigridSolver<Dim> msolver(bc, "FMG");

    int N = IntOp::power(n + 1, Dim);

    // 构造精确解
    Grid<Dim> x(n + 1, 0);
    do
    {
        x.value() = Real(u(domain.gridPoint(x)));
    } while (x++);
    // cout << x << endl;

    Tensor<Dim> v(n + 1, 0);
    msolver.solve(v);

    // 计算误差范数
    Real norm = TensorOp::grid(x - v, 2);
    cout << norm << endl;

    return norm;
}

void test()
{
    int k = 5;
    Vector x(k, 0);
    int n = 4;

    for (int i = 0; i < k; i++)
    {
        x[i] = testPoisson(n);
        n *= 2;
    }

    cout << endl;
    cout << "Rate:" << TensorOp::rate(x) << endl;
}

int main()
{
    test();
    // testPoisson(4);

    return 0;
}